Entropic origin of stress correlations in granular materials 
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We study the response of granular materials to external stress using experiment, simulation, 
and theory. We derive an entropic, Ginzburg-Landau functional that enforces mechanical stability 
and positivity of contact forces. In this framework, the elastic moduli depend only on the applied 
stress. A combination of this feature and the positivity constraint leads to stress correlations whose 
shape and magnitude are extremely sensitive to the applied stress. The predictions from the theory 
describe the stress correlations for both simulations and experiments semiquantitatively. 



A striking feature of dry granular materials and other 
athermal systems is that they form force chain networks 
in response to applied stress, such that large forces are 
distributed inhomogeneously into hnear chain-hke struc- 
tures P, 13] ■ A number of experimental studies have visu- 
alized and quantified these networks in granular systems 
using carbon paper Q and photoelastic techniques 0, Q . 
These studies demonstrated that geometrical and me- 
chanical properties of force chain networks are acutely 
sensitive to preparation procedures, especially near the 
jamming transition^. For example, in isotropically com- 
pressed systems, force networks are ramified with only 
short-ranged spatial correlations of the stress. In con- 
trast, in sheared systems, aligned force chains give rise to 
longer ranged stress correlations in the compressive direc- 
tion. Since granular (or other) systems near jamming are 
fragile and highly sensitive to preparation, one expects 
that their mechanical properties near jamming might not 
be captured by simple linear elastic response 

Developing alternative theoretical descriptions for 
granular media is challenging for several important rea- 
sons 8, 9, 10, JJJ: (1) since tensile stresses are absent 
in dry granular materials they only remain intact via 
applied stress, making the limiting zero-stress isostatic 
state, where the number of degrees of freedom matches 
the number of constraints [l3. 13 . [l^ . problematic; (2) 
forces at the microscopic level are indeterminant due to 
friction and disorder; (3) granular materials are athermal, 
so that conventional energy-based statistical approaches 
are not appropriate; and (4) near isostaticity we expect 
fluctuations to be important, both within a single realiza- 
tion of a system and from realization to realization. New 
methods are needed to bridge the gap between force net- 
works at small length scales and continuum elasto-plastic 
theory at large scales, and to capture the highly sensitive, 
fluctuating behavior of granular systems near jamming. 

We construct a model for stress fluctuations based on 
grain-scale force and torque balance and positivity of con- 
tact forces rather than energy conservation. We then 
calculate stress correlations and predict differences for 



systems under isotropic compression versus shear stress. 
We also perform complementary numerical simulations 
and experimental studies of jammed granular systems in 
2D subject to isotropic compression and pure shear. The 
stress correlation functions from theory, simulation, and 
experiment are in qualitative and in some cases quanti- 
tative agreement. In particular, the theory predicts that 
the form of the stress correlations depends on how the 
jammed states were prepared. 

Theoretical Framework: Fluctuations are inher- 
ently related to the number of microscopic states avail- 
able under a given set of macroscopic conditions. In equi- 
librium thermodynamics, the microcanonical entropy de- 
scribes the nature of fluctuations and response. When we 
turn to granular systems, the identification of states by 
energy is no longer useful, and a different criterion for 
identifying states is needed. The approach that we pur- 
sue here exploits a different conservation principal, based 
on force and torque balance, which applies rigorously for 



granular materials [1 



The force-moment tensor of mechanically stable (MS) 
packings, T, — J d'^ra{r), where a{r) denotes the local 
stress tensor, is an extensive variable that is a topolog- 
ical invariant [3, In the force- moment ensemble, 
S remains fixed barring system-spanning changes in a. 
Hence local fluctuations and response only involve grain 
configurations with the same E. To construct a theory for 
stress correlations, we adopt a coarse-grained approach, 
in which jammed configurations are represented by a con- 
tinuous field [2l[ and the entropy S{Y,) is defined via 
an appropriate Ginzburg-Landau functional. The theory 
should be valid close to the jamming transition where 
grains have negligible deformations, and stress fluctua- 
tions decouple from volume fluctuations. The decoupling 
is exact only in the limit of infinitely rigid grains [l8| . 

In both two and three dimensions, a continuous field 
can be defined that upholds force and torque balance 
[20I 1 of granular packings. We will focus on 2D systems, 
where a scalar field the Airy stress function [7], is 
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FIG. 1: Pressure correlations in Fourier space for pure com- 
pression. Angle-averaged S{q)K{r) from (a) simulations at 
different F (reduced units) and static friction coefficient ^, 
(b) experiments at different F (in units of N ■ m), /i — 0.7 
and z ranging from 3.35-3.68. The experimental data have 
been scaled by K{r) = {ziso/2 + c{z - zl^))/r'^^, with 
c — 2.8. Theory predicts that S{q)K{r) is independent of F 
for q < 1/C 



related to the local stress tensor by 



a(r) 



(1) 



We define T = TrE and r = Si - S2 = ^/T^ -4(detS)2, 
where si > S2 are the eigenvalues of E. Given a E, ^' 
can be expanded as ^'o + ipj where ip represents fluctu- 
ations around VPo- The field ^0 satisfies the biharmonic 
equation^-^^o = 0, and E^j = / d''r((5y V^^o - 

The probability for fluctuations ip can be written as 
= Z-i(E)e-'-[*o,^] ^ Z-i(S)e-'-fi['^l. The func- 
tional L|,[?/;] measures the contribution to the entropy 
from all grain packings that have a coarse-grained repre- 
sentation '4>{r). The partition function Z(E) = e"^*-^-* — 



generates correlators of the field ■0 [2]J . Be- 



cause of gauge freedom, can only depend on sec- 

ond derivatives of ip. Independent second order scalars 
involving second derivatives of ip can be constructed from 
the invariants of the local stress tensor, Tr(o-) and det{a). 
The coefficients of these terms and, therefore, the nature 
and strength of fluctuations are controlled by the field 
^'o (or E). To lowest order in -0, Lg['0] resembles the free 
energy for an elastic material in two dimensions 



+a3(t) {d.,dyi^f + ai{t) {dl4>){dl4>)] ■ (2) 

The crucial differences between the description in Eq. [2] 
and traditional elasticity theory are (a) the stiffness con- 
stants are determined by \E'o or E, and thus the theory 
is inherently nonlinear, and (b) the origin of L|,[7/;] is en- 
tropic. The functional \-f\'4'\ can be used to calculate av- 
erages and correlation functions. Below, we consider the 
behavior of pressure correlation functions under isotropic 
compression and pure shear using Eq. ((2). 



Isotropic Compression: For isotropic compression 
with no deviatoric stress, the stiffness constants only de- 
pend on r, and Lr['0] is isotropic: 



A(r)(9,a^V) 



{dim 



with two stiffness constants K and A. The positivity 
of contact forces implies that both d^"^ and d^"^ must 

This is a difficult constraint to 



be non-negative [15|, l2C 
impose exactly on the fluctuating fleld ■0. However, it 
can be enforced in a mean-fleld way by requiring that 
K{T) > l/r^, which guarantees that the amplitude of 
the long- wavelength fluctuations are such that the posi- 
tivity criterion is met 2^, 2^- The inequality constraint 
on K{T) implies that different preparation histories can 
lead to different fluctuations. The maximum entropy of 
jammed packings is achieved in protocols that meet the 
equality, and we focus here on these marginal packings. 
Note that the positivity constraint does not impose any 
conditions on A, which is therefore taken to be indepen- 
dent of r. Near jamming when F ^ 0, K(T) becomes 
arbitrarily large and the A terms can be ignored . 

The results for the correlations of the local pressure 
are best visuahzed in Fourier space. From Eq. ([3]), these 
correlations are predicted to be isotropic: 



5(q) = <|(5r(q)|2>=q4<|VXq)P> 

1 + eq^ ' 



(4) 



where ST = V^-;/;, q is the wavevector, and ^ is a cor- 
relation length that describes the decay of correlations 
at large q, and is defined by higher order terms not in- 
cluded in Eq. ([3]). In an experiment or simulation at fixed 
T/A, there are many MS packings, and each is charac- 
terized by a continuously varying field ip(j)- The spa- 
tial correlations of stress, for a given F, are determined 
by averaging over these configurations. If the configura- 
tions are sampled according to the theoretically predicted 
P[tp] oc e~'"'"['^l, the correlations measured in simulations 
and experiments should be well described by the field- 
theoretic predictions. Since frictionless granular pack- 
ings are isostatic near jamming t24|, an exact calcula- 
tion yields K{T) = Ziso/(2F2) [20|, where the number 
of contacts Ziso = 4 in 2D. In contrast, frictional pack- 
ings have Ziso = 3 in 2D, but are often hypostatic [2^. 
We do not have an exact result for K{r) away from iso- 
staticity, however, a form that has been successful [20|] is 
K(r) — (ziso/2 + c(z — z?q))/F2, where c is a phenomeno- 
logical constant. We will use this form to compare the 
predictions with results from frictional packings. 

To test the g-space pressure fluctuations predicted in 
Eq. we have numerically generated MS packings 
of bidisperse disks {N/2 large and N/2 small particles 
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with diameter ratio r = 1.4) both with and without 
friction near the jamming transition using well-known 



27| . For frictionless 



packing-generation algorithms [26 
grains, these algorithms generate packings at the margin 
of stability [1^]. We studied system sizes ranging from 
N = 256 to 4096, systems with square cells and periodic 
boundary conditions, pressures in the range T/A — 10~^ 
to 10~^ (in reduced units of the grain stiffness), and static 
friction coefRcients in the range /i = [0, 1]. We have also 
carried out experiments using a biaxial apparatus that 
has been described previously [2^, [2§|. The biax is 
a device that allows us to apply highly controlled defor- 
mations to quasi-2D systems of photoelastic disks. By 
using photoelastic disks, it is possible to obtain all con- 
tacts and contact forces in the system. In this study, con- 
tact forces are calculated for N — 1228 disks, with 
large and AN/5 small disks with diameter ratio r — 1.2 
and coefficient of static friction /i = 0.7. The experimen- 
tal protocol generates packings farther from isostaticity 
than those from the simulation protocol. 

Results for pressure correlations in compressed systems 
are shown in Fig. [TJ In general, we find that S{q) decays 
isotropically with q. In Fig. [TJa), we plot the angle- 
averaged S{q) from simulations, normalized by K{r) at 
r/A = 10~^ and 10""* as a function of qD, where D is the 
small particle diameter. For both frictional = 1 and 
frictionless /i = grains, the results from the simulations 
match Eq. (|4]) at small q, with no fitting parameters. 
In Fig. [TJb) we plot the angle-averaged S{q)K{T) from 
experiments. Both simulations and experiments confirm 
the theoretically predicted scaling of S{q) as 1/K{r) and 
increasing stiffness as the system is decompressed. 

Pure Shear: In the presence of an imposed pure 
shear, the positivity constraints on the stress lead to dif- 
ferent conditions on the stiffness constant K in the x and 
y directions. A dramatic consequence is that the pressure 
correlations S{q) become anisotropic even for infinitesi- 
mal shear, and the correlations in real space become long- 
ranged. To lowest order in ip the entropy functional is 



\(d,dyi,f) + K'{T,T)dl4>dli^}. 



(5) 



Here, x (y) is the principal axis of E with the smaller 
(larger) eigenvalue. In the case of pure shear, there are 
now two distinct stiffness coefficients to ensure that both 
axx and ayy are positive. In addition, K' controls the 
entropy cost of fluctuations that contribute to Oxx^yy A 
more symmetric version can be constructed by demand- 
ing that K' ^ yJ{K{T + t)K{T -t)), although we have 
no rigorous argument to support this form. 

The pressure correlations predicted from Eq. ([5]) are: 



5(q) = qHK{T + T)qt, + K{T-T)q^y + 2K'{T)qlql 

+ K{v)eq^)-\ 
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FIG. 2: Contours of S{q)K{T) under pure shear: (a) theory, 
using T /T = ^/D = 0.3; (b) experiment, with intensity on a 
log-scale; and simulations of (c) frictionless and (d) frictional 
(/I — 1) particles. In both sets of simulations, r/F — 0.3. In 
all plots, compression (dilation) is along the vertical (horizon- 
tal) axis. 



The anisotropic, dipolar nature of this correlation func- 
tion is depicted in Fig. [2] (a). To compare theory with 
experiment, we create a sheared packing by first isotrop- 
ically compressing the system to a mechanically stable 
state at a density slightly above jamming. We then ap- 
ply pure shear by expanding the system in one direction 
while compressing in the other, keeping the density con- 
stant. The resultant pressure correlations are shown in 
Fig [21 (b), and they match the expected form within the 
noise of the data. To compare theory and simulation, we 
generated MS packings of bidisperse disks with and with- 
out friction over a range of stress ratios r/F. To do this, 
we compressed (dilated) the simulation cell in the y (x) 
direction by e = 6L/L over the range e — [10^^, 10^'^]. 
Pressure correlations from simulations in Fig. [2] (c) and 
(d) also show the dipolar character of the pressure cor- 
relations. A key prediction of Eq. [Slis that limg^o 5'(q) 
depends on the direction of approach. This feature is 
clearly demonstrated in Fig. [3l which shows the simu- 
lation and experimental results for iS'(q) along different 
cuts in g— space, along with the small-g predictions from 
theory. There is good agreement between theory and 
simulation for q^ — and = 0, where theoretical pre- 
dictions exist. Even though the theoretical predictions 
make several simplifying assumptions such as z — Zieo 1 
(small F) and r/F <C 1, we observe qualitative agreement 
with experimental data. In particular the pressure cor- 
relations depend on the direction of approach to g = 
and they are larger along qy = than q^ — 0. 

The anisotropic nature of correlations in g-space imply 
anisotropic decays in real space [SQ] with a slower decay 
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FIG. 3; Cuts along axes specified in the legends for S{q)K{r) 
under pure shear, (a) Simulation results with fi — and 
r/r = 0.3. The solid lines are theoretical predictions for 
qx ~ and = 0. Results for jj, — 1 are qualitatively 
similar, (b) Experimentally measured S{q)K{T) contours at 
r/r = 0.51. The ratio of S{qy = 0)/S{qx = 0) is close to the 
theoretical prediction in (Eq. |6j, 1 + 4^ + 0((f )^) ~ 3. 



along the direction of higher compression. The entropic 
formulation with the positivity constraint, therefore, pro- 
vides an explanation for the shear-induced anisotropy in 
pressure correlations observed in experiments f^. 

Discussion: We present a field theoretic approach, 
based on entropy of packings, for describing stress fluc- 
tuations in granular packings. The theory enforces con- 
ditions of mechanical stability and positivity of contact 
forces, and applies close to the jamming transition, where 
grains have small deformations. From the theory, we cal- 
culate pressure correlations and show that they depend 
sensitively on the method used to generate the jammed 
state. Under isotropic compression, all correlations are 
isotropic and obey a simple scaling relation as a function 
of compression. For packings subjected to pure shear, the 
correlations are anisotropic with a characteristic dipolar 
feature in g-space. The anisotropy is a consequence of the 
positivity constraint, which causes g-space stress fluctua- 
tions to be reduced along the compressive direction. The 
present approach provides a means of relating stress fluc- 
tuations to the history of granular systems, which deter- 
mines the force moment tensor, and an explanation for 
the anisotropic behavior of stress fluctuations. The the- 
oretical predictions for the pressure correlation functions 
are confirmed, semiquantitatively, by simulations of MS 
packings with and without friction and by experiments 
on photoelastic disks. This agreement is remarkable since 
it validates the idea that the entropy of MS packings can 
be used to determine the response of the this far-from- 
equilibrium system. 

Work supported by "nsf-dmr0555431 (BC,SH), 
bsf-dmrn4 48838 (GL), nsf-dms0835742 (CO), 
nsf-dmr0555431 (JZ,TSM,RB), and YINQE (GL). 
BC acknowledges discussions with Nick Read, and CS, 
BC, GL acknowledge the Aspen Center for Physics 
and Lorentz Center, where aspects of this work were 
performed. 



[2: 

[s: 
[4: 

[6; 

[7: 

[8 

[9: 

[10: 
[11 

[12 

[13 
[14 

[15 

[16 

[17 
[18 
[19 
[20 
[21 

[22' 
[23 

[24 

[25 

[26 

[27 

[28; 

[29; 

[3o; 



C. -h. Liu, S. R. Nagel, D. A. Schecter, S. N. Coppersmith, 
S. Majumdar, O. Narayan, and T. A. Witten, Science 
269, 513 (1995). 

J. Zhou, S. Long, Q. Wang and A. D. Dinsmore, Science 
312, 1631 (2006). 

D. M. Mueth, H. M. Jaeger, and S. R. Nagel, Phys. Rev. 
E 57, 3164 (1998). 

P. Dantu, Ann. Ponts Chaussees 4, 144 (1967). 
J. Ceng, D. Howell, E. Longhi and R. P. Behringer, G. 
Reydellet, L. Vanel, E. Clement, and S. Luding, Phys. 
Rev. Lett. 87, 035506 (2001). 

T. S. Majmudar and R. P. Behringer, Nature 435, 1079 
(2005). 

L. D. Landau and E. M. Lifshitz, Theory of Elasticity, 

(Butterworth-Heinemann, London, 1986). 

M. Otto, J.-P. Bouchaud, P. Claudin, and J. E. S. Soco- 

lar, Phys. Rev. E 67, 031302 (2003). 

M. E. Gates, J. P. Wittmer, J.-P. Bouchaud and P. 

Claudin, Phys. Rev. Lett. 81, 1841 (1998). 

R. Blumenfeld, Phys. Rev. Lett. 93, 108301 (2004). 

I. Goldhirsch and C. Goldenberg, Eur. Phys. J. E 9, 245- 

251 (2002). 

A. V. Tkachenko and T. A. Witten, Phys. Rev. E 60, 
687 (1999). 

C. F. Moukarzel, Phys. Rev. Lett. 81, 1634 (1998). 

M. Wyart, S. R. Nagel, and T. A. Witten, Europhys. 

Lett. 72, 486 (2005). 

S. Menkes, C. S. O'Hern and B. Chakraborty, Phys. Rev. 
Lett. 99, 038002 (2007). 

R. Blumenfeld in Lecture Notes in Complex Systems Vol 
8: Granular and Complex Materials, edited by T. Aste, 

A. Tordessilas and T. D. Matteo (2007). 

B. P. Tighe, A. R. T. van Eerd and T. J. H. Vlugt, Phys. 
Rev. Lett. 100, 238001 (2008). 

C. Song, P. Wang and H. A. Makse, Nature 453, 629 
(2008). 

R. C. Ball and R. Blumenfeld, Phys. Rev. Lett. 88, 
115505 (2002). 

S. Henkes, Ph.D. dissertation, (2008); S. Henkes and B. 

Chakraborty, to appear in Phys. Rev. E. (2009). 

N. Goldenfeld, "Lectures on Phase Transitions and the 

Renormalization Group", (Addison- Wesley, New York, 

1992). 

B. Tighe, private communication. 

Higher-order derivatives of will likely enter if A' ^ 00 
and suppress stress fluctuations. 

G. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, 
Phys. Rev. E 68, 011306 (2003). 

L. E. Silbert, D. Ertas, G. S. Crest, T. C. Halsey and D. 
Levine, Phys. Rev. E 65, 031304 (2002). 
N. Xu, J. Blawzdziewicz and C. S. O'Hern, Phys. Rev. E 
71, 061306 (2005). 

H. P. Zhang and H. A. Makse, Phys. Rev. E 72, 011301 
(2005). 

T. S. Majmudar, M. Sperl, S. Luding and R. P. 

Behringer, Phys. Rev. Lett. 98, 058001 (2007). 

J. Zhang, T. S. Majmudar, A. Tordesillas, and R. P. 

Behringer, preprint (2009). 

S. Henkes, unpublished (2009). 



